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Many-body QCD in leading high energy Regge asymptotics is described by the Balitsky-JIMWLK 
hierarchy of renormalization group equations for the x evolution of multi-point Wilson line corre- 
lators. These correlators are universal and ubiquitous in final states in deeply inelastic scattering 
and hadronic collisions. For instance, recently measured di-hadron correlations at forward rapidity 
in deuteron-gold collisions at the Relativistic Heavy Ion Collider (RHIC) are sensitive to four and 
six point correlators of Wilson lines in the small x color fields of the dense nuclear target. We 
evaluate these correlators numerically by solving the functional Langevin equation that describes 
the Balitsky-JIMWLK hierarchy. We compare the results to mean- field Gaussian and large iVc ap- 
proximations used in previous phenomenological studies. We comment on the implications of our 
results for quantitative studies of multi-gluon final states in high energy QCD. 
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I. INTRODUCTION 

QCD in high energy Regge asymptotics can be de- 
scribed as a dense many-body system of "wee" gluons 
and sea quarks. In the infinite momentum frame, glu- 
ons with transverse momenta k± < Qs saturate phase 
space maximally, where Qs{x) is a dynamical saturation 
scale [1] that grows with decreasing fractions x of the lon- 
gitudinal momentum of the hadron carried by the gluons. 
The properties of saturated gluons are described by the 
Color Glass Condensate (CGC) effective theory [21 |3], 
where the degrees of freedom are static color sources in 
the hadron at large x, coupled to the dynamical wee gluon 
fields at small x. Renormalization group equations, de- 
rived from requiring that observables be independent of 
the separation in x between sources and fields, lead to an 
infinite hierarchy of evolution equations in x, for n-point 
Wilson line correlators averaged over dense color fields in 
the hadron. Given appropriate initial conditions at large 
X, solutions of this Balitsky-JIMWLK hierarchy [4", al- 
low one to compute a wide range of multi-particle final 
states in deeply inelastic scattering (DIS) and hadronic 
collisions. 

A prominent example is provided by inclusive DIS 
structure functions F2 and Fl^ which are proportional 
to the forward scattering amplitude of a qq "dipole" on 
a nucleus. The forward dipole amplitude (dipole cross 
section) can be expressed as 

;i-iTvr(b.+i)r.(b.-f)). 

where = xt — yr is the transverse size of the dipole, 
bT = (xT+yT)/2 is the impact parameter relative to the 



hadron, and the rapidity Y = ln(xo/x), where xq is the 
initial scale for small x evolution. The dipole amplitude 
is the expectation value D = {D) of the dipole operator 
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I)(XT -yT) = ^ Tr F(xT)Ft(yT) 



(2) 



averaged over the color sources of the target evaluated 
at the rapidity Y . This average obeys the Balitsky- 
JIMWLK equation that relates its energy dependence to 
the expectation value of a four-point operator: 



d . . TVcQ^s 

D{^T - Yt) 



dY 
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(xt - yr)^ 



(xt - zt)^(zt - Yt^ 



X (D{yiT - zt) D{zt - Yt) - D{yiT - Yt)) ' (3) 



In the large Nc approximation the expectation value of 
Z)^ factorizes and the equation becomes a closed one. 
This is known as the Balitsky-Kovchegov (BK) equa- 
tion HIS]: 



:D[^T - Yt) = 



[ (xt - Yt)^ 

dF" "'^ 2^2 J^^ _ zt)^zt - Yt)^ 

X [L)(xt - zt) D{zt - Yt) - D{^t - Yt)] • (4) 



In addition to Eq. ([T]), this dipole correlator appears 
in a number of final states in both DIS and hadronic 
scattering; the BK equation for its energy evolution is 
widely used in phenomenological applications. The mean 
field approximation (F)'^) ~ (D)'^ has been checked by 
numerical solutions of the JIMWLK equations, and it 
is seen that it is satisfied to a very good approximation 
(much better than the one might expect) [ZlIH]. 

Unless noted otherwise, throughout the paper we define 



2 



the saturation scale Qs(^) from the expectation value of 
the dipole operator as D (r = \/2/Qs) = e~^/^. 

For less inclusive observables, new universal degrees 
of freedom beyond dipoles are encountered. Examples 
include small-x di-jet production in e+A DIS quark- 
antiquark pair production in hadronic collisions [10] and 
near-side long-range rapidity correlations [TT. Here we 
focus on n-point functions which appear in forward di- 
hadron production in light on heavy hadron collisions, 



p -\- A — > hi h2 X, di process that has been studied re- 
cently in deuteron-gold collisions at RHIC. When both 
hadrons are produced at forward rapidities in the pro- 
ton/deuteron fragmentation region, the dominant under- 
lying QCD process is the scattering of a large Xi va- 
lence quark from the deuteron off small X2 partons in 
the nuclear target, with the emission of a gluon from the 
valence quark either before or after the collision. This 
cross-section is expressed as [9l [12] , 



OC 

d3fci d3fc2 2 



^XT,XT,yT,yT \ 



Z)(yT,XT)^(xT,ZT)-^(zT,XT)^(xT,yT) + ^^(zT,ZT) + ^ (^Z)(yT, Zt) +^(zt, Yt) -^(YT, Yt)) \ , (5) 



with zt = z^T + (1 — ^)Yt and likewise, z^ = z5tT + 
(1 — 2;)yT- ^ denotes the splitting function for producing 
a photon off a quark, with the four co-ordinates denot- 
ing the transverse spatial co-ordinates of the quark and 
gluon in the amplitude and conjugate amplitude. The 
color field dynamics specific to gluon emission are ab- 
sorbed in the expectation value (), which contains a new 
quadrupole operator. 



(5(xt, Yt,ut, vt) 



TrF(xT)Ft(yT)F(uT)Ft(vT). 

(6) 

We denote the expectation value of the quadrupole op- 
erator by Q = (Q). Unlike the (Z)^) in Eq. (|3|, it is not 
reducible to the product of dipoles even in the large 
and large A approximations and is a novel universal cor- 
relator in high energy QCD 13 , interesting both from 
theoretical and phenomenological perspectives. 

In this paper, we determine the expectation values 
of relevant Wilson line correlators for a SU(3) gauge 
group explicitly numerically. It is known that the evolu- 
tion of the expectation values of Wilson line correlators 
can be expressed as a functional Fokker-Planck equa- 
tion [Hj, which in turn can be re-expressed as a func- 
tional Boltzmann-Langevin equation for the Wilson lines 
themselves [15], 



dF(xT) 

dr 

where 



(7) 



e{^T,ZT), =[-) [1-U{^tVU{zt)] 

(8) 

is the "square root" of the JIMWLK Hamiltonian. Here 
t/'s are adjoint Wilson lines, which are related to fun- 
damental Wilson lines through the identity, Uabi^r) = 



2Tr(ri/t(x^)tV(xT)). The equation includes a term 
proportional to a Gaussian white noise ^ satisfying 
(C(xt)?) = and (C(xT)??(yT)5) = S^'S'^S^^H^t - Yt)- 
Finally, there is also a drag term 



<t(xt)" 
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(xt - Zt)2 



Tr [r''[/(xT)tC/(zT)] 



(9) 



where is a generator of the adjoint representation. 



II. NUMERICAL METHOD 

As alluded to previously, numerical solutions of the 
Boltzmann-Langevin hierarchy were obtained for fixed 
coupling [7 and used to study the factorization of dipole 
correlators. Apart from the running coupling, we use 
the numerical method of ref. [7 for solving the JIMWLK 
equation. Several running coupling prescriptions have 
been suggested in the literature [16]; in this paper, as in 
[17j . we will assume that the coupling constant runs as 
a function of the "daughter" dipole size r = |xt — zt|. 
The Landau pole is regulated by taking 



47r 



(10) 



with a parameter c which regulates the sharpness of the 
cutoff. The scale A in the coupling is parametrically of 
the order of Aqcd^ but the exact value that should be 
used is scheme dependent. 

The initial conditions are those of the McLerran- 
Venugopalan (MV) model [18 ; one has for the initial 

) with 



rapidity Y = 0, V{xt) = Ofli exp( 
{pt{^T)p1{yT)) 



: ^/V^^^^H^T - yT)S''^S^^/Ny, with the 
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indices k^l = 1 . . . Ny representing a discretized longitu- 
dinal coordinate, taking care of the finite extension of the 
source in the x~ direction. The normalization is chosen 
such that Efc,,{Pfe(xT)p?(yr)> = 9^^l^S(^^^T - y^)<5°^ 
Given the numerical implementation of this initial con- 
dition [19] for the V's, using a longitudinal resolution 
of Ny = 100, one can then solve Eq. ([7| on a 2-D lat- 
tice. The Poisson equation is solved by leaving out the 
zero transverse momentum mode. This procedure cor- 
responds to an infrared cutoff given by the size of the 
system. The calculation is performed on a regular square 
lattice of sites with periodic boundary conditions. Its 
volume is {Nto)'^ where a is the lattice spacing. In the 
calculations presented we use Nt = 512 unless otherwise 
stated, and 8s = 0.00026, g'^/ia = 0.109375. The pa- 
rameters controlling the running coupling are taken as 
c = 0.2, A = 0.0536 ^^/i, and /io = 2.5A. These param- 
eters are chosen to be close to the phenomenologically 
realistic range for the speed of evolution as observed in 
fits to F2(x, Q^) data, but have not been adjusted to give 
a best possible fit. 



On inspection, it is apparent that this approx- 
imation is problematic because it does not re- 
duce to all the right "coincidence limits"; tak- 
ing ut = vt one has Q(xt, yr, ut, ut) = 
Dip^TiyT)-) but the r.h.s. of Eq. ( pTj ) reduces to 
I (I)(xt, Yt) + ^(xt, ut)^(ut, Yt)) instead. Plots of 
the JIMWLK solution for the line and square configura- 
tions respectively compared to the "naive" approxima- 
tion 



are shown in Fig. [2] It is apparent that the approxima- 
tion fails even for the initial condition at F = 0, a result 
that is not ameliorated by the JIMWLK RG evolution. 
The disagreement between the JIMWLK and approxi- 
mate results is greater for the line configuration than the 
square configuration. 



III. NAIVE LARGE APPROXIMATION 
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xr 
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FIG. 1: Square arrangement of the coordinates in Eq. ([6|. 
The filled circles represent Wilson lines and the open ones 
conjugates. 

We will first present results for the expectation value 
of the quadrupole correlator Q. Since this correlator is 
a function of four independent two dimensional spatial 
vectors, we will for simplicity study its properties for two 
specific spatial configurations, one of which is the square 
configuration Qu(j)i which has the four coordinates ar- 
ranged in a square of size r = |xt — yrl as shown in 
Fig. [ij The other is a simple line configuration Q\(r)^ 
where ut = xt, vt = yr and r = |xt — yrl = |ut — vt|. 
The expectation values of all correlators are computed by 
averaging over all positions on the 2D lattice and 20 dif- 
ferent initial configurations. 

A "naive" large approximation for Eq. (|6| consid- 
ered previously in phenomenological studies is 



IV. GAUSSIAN APPROXIMATION 



The fact that the expectation value of the quadrupole 
correlator does not factorize in the naive way of Eq. (11) 
was pointed out in [20 and is seen explicitly already in 
the MY model that specifies the initial conditions^. We 
will in this paper compare the JIMWLK result for Q 
to a "Gaussian" approximation [22l |23] (also referred 
to as Gaussian Truncation [3[ [8]). This is obtained 
by assuming that the correlators of the color charges 
are Gaussian variables even after JIMWLK evolution, 
and therefore all the higher point functions can be ex- 
pressed in terms of a single two point correlation func- 
tion, — ln(I)(xT — yr)) = ^r(xT — yr) in the notations 
of ref. [9 . We obtain this two point function from our 
solution of the JIMWLK equation. However, as shown 
in ref. [8 , using the solution of the BK equation for the 
two point function would also be a very good approxima- 
tion. The Gaussian approximation has also been moti- 
vated formally at asymptotically small x in ref. [24] . Thus 
in the Gaussian approximation the higher point functions 
are related to the two point function similarly as in the 
MY model. For the specific case of the quadrupole Q the 
explicit expression has been derived in ref. [9]. 



For the square and line configurations we consider 
here, the cumbersome general Gaussian expression for 



Q(xT,yT,UT, Vt) 



1 



(i:>(xT,yT)^(uT, Vt) 

- i^(xT, vt)1^(ut, Yt)) (11) 



^ We leave the study of quartic contributions [2j to the MV action 
for future study. 
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FIG. 2: The four point function Q for coordinates in the Une configuration (left) and the square configuration (right) for the 
initial condition {Y — 0) and after evolution {Y — 5.2). The JIMWLK quadrupole expectation value (solid lines) is compared 
to the naive approximation in Eq. (12) (dashed lines). 
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FIG. 3: The four point function Q for coordinates in the line configuration (left) and the square configuration (right) for the 
initial condition {Y — 0) and after evolution {Y = 5.2). The JIMWLK quadrupole expectation value (solid lines) is compared 
to the Gaussian approximation (dashed lines) Eqs. (13), (14) . 



the quadrupole correlator simplifies greatly and we find, 



Q\{r) 



Nr. 



■ Nc + 1 



-(D{r)Y 

^(^MJ 



Qn(r) « {D{r)) 

iVc<00 



A^'c + 1 / D{r) \ 



(13) 



(14) 



2 \D{^/2r) 

iVc-1 / D(v^r) \^- 
2 1^ D{r) ) 

In Fig. |3j the numerical results from the solution of 
the JIMWLK RG equation for the quadrupole (in the line 
and square configurations) are compared to this Gaussian 
approximation. When computing the Gaussian approxi- 
mation, we chose to only calculate the averages in D from 
contributions aligned with the square or the line config- 
uration, respectively. While the agreement of Q with its 
Gaussian approximation for the initial condition F = is 



required by definition (because one has MV initial con- 
ditions in both cases), the agreement of the JIMWLK 
simulation with the evolved Gaussian approximation is 
remarkably good. This suggests that the computations 
in ref. [22| ^25^ that rely on this approximation may, at 
least in this aspect, be robust. 

If one takes the large A^c limit of the expressions in 
Eqs. (flsl) and (IT4|), one finds that 



Q|(r) = D\r)[1^2\n{D{r))] 

Nc^oo 



Qair) ^= D\r) 



1 + 21n 



D{y/2r) 



(15) 
• (16) 



The JIMWLK result compared to this particular large 
A/'c approximation is shown in Fig. [4j The agreement of 
the large Gaussian approximation with the JIMWLK 
result is quite good for Qn, especially for the initial condi- 
tion; it is less so for Q\. In the latter case, the agreement 
for small rQs is good, with discrepancies showing up for 
rQs > 1. We note that the difference between the "naive" 
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FIG. 4: The four point function Q for coordinates in the Une configuration (left) and the square configuration (right) for the 
initial condition {Y — 0) and after evolution {Y — 5.2). T he J I MWL K quadrupole expectation value (solid lines) is compared 
to the large A^c limit of the Gaussian approximation, Eqs. (15), (16) (dashed lines). 



approximation ^12]) and this Gaussian large A^c approx- 
imation (13), (|14[) is in the additional logarithmic term 



in the latter. In the case of the line configuration, this 
provides some insight into the discrepancy of the naive 
approximation with the result from JIMWLK evolution. 
The latter is not constrained to be greater than zero, 
while the former is. The additional logarithmic term re- 
laxes this constraint because it can change the sign of the 
result. 



the dipole and A^^,A^'^ for the quadrupole in the line 
and square configurations respectively. The plot shows 
that initially the quadrupole evolves more rapidly than 
the dipole. After about 6 units in rapidity, the evolution 
reaches a universal geometrical scaling regime and the A 
parameters settle to a common value. 



VI. SIX POINT FUNCTION AND DIHADRON 
CORRELATIONS 



V. GEOMETRIC SCALING OF THE 
QUADRUPOLE 

We now turn from this comparison of the quadrupole 
expectation value to study aspects of its evolution deter- 
mined from the solution of the JIMWLK equation. The 
KG evolution of the quadrupole has been studied pre- 
viously ^31 [26] in the large limit of Mueller's dipole 
model ^27]. In particular, it has been argued very recently 
in ref. [26] that quadrupole evolution should demonstrate 
geometrical scaling [28 . 

Our results are shown in Fig. [5j On the left we plot 
Qn as a function of rQs- We observe that after initial 
transient behavior the amplitude Qn settles on a shape 
that is a universal function of rQs- Thus for a given 
X and probed in a process, the quadrupole ampli- 
tude depends only on the combination proportional to 
Qs{xY IQ'^ thereby demonstrating geometrical scaling for 
this quantity. 

We previously defined the usual saturation scale 
Qs through the dipole operator as D (r = a/2/Qs) = 
e~^/^. One can analogously characterize the evolution 
of the quadrupole by introducing "quadrupole saturation 
scales" corresponding to the square and line coordinate 

arrangements. We define these as (5|,n = a/2/Qs'^^ = 

e~^/^. On the right of Fig.[5]we show the evolution speeds 
A = dlnQs/ dY of these different saturation scales, A for 



Let us now return to the expression we had in Eq. ([5| 
for forward di-hadron production in hadronic collisions. 
Experiments at RHIC for deuteron-gold scattering at 
high energies have shown that the away-side peak in 
di-hadron correlations is significantly broadened for cen- 
tral collisions at forward rapidities [29 as predicted in 
the CGC framework [12] and confirmed in more detailed 
analyses [30 . However, these analyses relied on factoriza- 
tion assumptions that, as we have seen, are not justified 
because the quadrupole correlator is not simply factor- 
izable. It is not the quadrupole correlator that appears 
directly in Eq. ([5|, so we shall now consider the expres- 
sion 



6'6(xT,yT,UT, Vt) 



N' - 1 



X (^Q(xT,yT,UT, vt)1^(vt,ut) - D{:siT .Jt) J , 

(17) 

corresponding to the first and last terms in the brack- 
ets in Eq. ([5|. The JIMWLK evolution equation for this 
quantity was derived recently [20 . As for the quadrupole 
correlator, we would like to compare numerical results 
for Sq to a Gaussian approximation. The latter, how- 
ever, would strictly require that one compute the prod- 
uct (QD) in this approximation. Since results for this 
quantity are not at present available, we will assume that 
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FIG. 5: Left: The JIMWLK quadrupole amplitude Qn for the square configuration versus the scaling variable QsV (on a 
logarithmic axis). After initial transient behavior, the amplitude settles to a universal curve (for Y > 7.8) which depends on 
QsV alone. Right: The evolution speed A = dlnQg/ dY extracted from the dipole amplitude D, the quadrupole amplitude in 
the two spatial configurations (line and square, A^^, A^^) and the six point function 5*6, also in the two spatial configurations 
(A^^,A^^). 
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FIG. 6: Left: The JIMWLK result for 5*6 (see text) plotted for the line configuration compared to the Gaussian approximation 
for Sg- Right: The JIMWLK result compared to a naive large A^c result. 



{QD) ~ {Q){D) and compare the Nc = 3 Gaussian ap- 
proximation for {Q){p) to the numerical JIMWLK re- 
sults for Sq. Figure [g] (left) shows the result for the line 
configuration for Sq. We observe that the agreement of 
the JIMWLK and approximate results is quite good for 
rQs ^ 1 or rQs ^ 1, but that there are noticeable de- 
viations in the region 1 < rQs ^ 3. Due to the good 
agreement of Q with the Gaussian approximation from 
Eqs. (13), ([T4| shown above, we interpret these devia- 
tions for Sg as 0{1/Nc) corrections to (QD) ^ {Q){D)^ 
For the square configuration of Sq (not shown) the devi- 
ations are much less, as was the case for the quadrupole. 
In Fig. [g] (right), we plot the JIMWLK results against the 
"naive" large Nc approximation 5'6|(^) ^ Sqc\{'^) ^ D{r)^ 
that has previously been considered in the literature. 
Once again the deviations are large, suggesting that this 
approximation is not tenable. We characterize the six- 
point function by the saturation scales Qs''^, defined by 



56,n(r = V^/Q^°) = e 



-1/2 



The corresponding evolu- 



tion speeds are also shown in Fig. [5] 

As a final result, we present a visualization of JIMWLK 
evolution that demonstrates the role of fluctuations. In 
high energy QCD, fluctuations from event-to-event can 
occur because of fluctuations in the impact parameter po- 
sitions of gluons, in their position in rapidity, and in the 
fluctuations in the number of gluons [31 . All of these 
fluctuations are captured in the numerical simulations 
of the JIMWLK hierarchy. These results are shown in 
Fig. [7| which shows the fluctuations of the Wilson lines 
in the transverse plane at different rapidities. The de- 
creasing of the correlation length ^ l/Q^ with energy is 
clearly visible. 

To summarize, we have in this work performed simula- 
tions of the running coupling SU(3) JIMWLK equation 
that describes the behavior of expectation values of Wil- 
son line correlators in high energy QCD. We have pre- 
sented first results for the evolution of specific higher n- 
point functions which are related to experimental observ- 
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FIG. 7: (animated online, requires Acrobat reader) Correlation 1 / Nc {V^ {0 , 0)V {x , y)) between the center position (0,0) and 
the point (x, y) for three different rapidities Y. This illustrates the degree of fluctuations and shows how the correlation length 
decreases dynamically with increasing Y. The first image can be animated to show the evolution with rapidity. 



ables in DIS and in p+p, p+A collisions. In particular, 
we find that the "quadrupole" correlator can be approxi- 
mated quite well by a careful Gaussian approximation for 
A^c = 3. The Nc ^ oo Gaussian approximation is accu- 
rate at short distances, rQg < 1, but may display more 
significant relative deviations in the saturation regime, 
rQs ^ 1- We have also provided evidence for travelling 
wave solutions and geometric scaling for the quadrupole. 
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